#### Does mislabeling COVID-19 elicit the perception of threat and reduce blame?
#### Chengxin Xu & Yixin Liu
#### JBPA Replication file: main tables and figures

# install.packages(c('list','ggplot2','plotly','ggpubr','egg','simpleboot','cowplot','Rmisc','list','xtable','stargazer'))

library(list)
library(ggplot2)
library(plotly)
library(ggpubr)
library(egg)
library(simpleboot)
library(cowplot)
library(Rmisc)
library(list)
library(xtable)
library(stargazer)

##import data
setwd()
exp<-read.csv("JBPA_Xu&Liu.csv")

####Subgroup Datasets####
exp.con<-subset(exp,con==1) #subgroup: conservatives
exp.lib<-subset(exp,lib==1) #subgroup: liberals
exp.mod<-subset(exp,mod==1) #subgroup: moderates
exp.dem<-subset(exp,Dem==1) #subgroup: Democrats
exp.rep<-subset(exp,Rep==1) #subgroup: Republicans
exp.ind<-subset(exp,Ind==1) #subgroup: Independents
exp.hs<-subset(exp,hsocial==1) #subgroup: High media use
exp.ls<-subset(exp,hsocial==0) #subgroup: low media use

####Figure Aesthetic####
myShape <- c(19,15,17)
pos2 <- position_dodge(width=0.2)
pos3 <- position_dodge(width=0.5)
limits <- aes(ymax = upr, ymin = lwr)

####Figure 2####
##List Item Count and Proportion 
##of Participants Who Reported Perceived Threat

# Figure 2 upper panel (item count)
ict.df<-summarySE(exp, measurevar="ict", 
                  groupvars=c("Item","ch1"),
                  conf.interval = 0.90)

# plot 2.1
f2.1.plot <- ggplot(ict.df, 
                    aes(y=ict, x=ch1, group=Item, shape=Item))+
  scale_shape_manual(values =c(1,16))+
  scale_x_discrete(labels=c("COVID-19","Chinese Virus"))+
  theme_article()
f2.1.plot <- f2.1.plot + geom_point(size=4, position=pos2) + 
  geom_errorbar(aes(ymin=ict-ci, ymax=ict+ci), width=0, colour="black", 
                position = pos2) + 
  labs(x=element_blank(),y="List Item Count") + 
  theme(legend.title=element_blank(),
        axis.title.y=element_text(size=10),axis.text=element_text(size=10),
        legend.justification=c(1,1), legend.position=c(1,1))+
  guides(colour=FALSE, shape = guide_legend(reverse=T))
f2.1.plot

# Figure 2 lower panel (proportion estimation)
# List estimation
st.mod.c = t.test(exp$ict[exp$group=="covid4"],
                  exp$ict[exp$group=="covid3"],
                  conf.level=0.90)
stc.est<-st.mod.c$estimate[1]-st.mod.c$estimate[2]
upper.stc<-st.mod.c$conf.int[2]
lower.stc<-st.mod.c$conf.int[1]
stc.ci<-upper.stc-stc.est

st.mod.t = t.test(exp$ict[exp$group=="virus4"],
                  exp$ict[exp$group=="virus3"],
                  conf.level=0.90)
stt.est<-st.mod.t$estimate[1]-st.mod.t$estimate[2]
upper.stt<-st.mod.t$conf.int[2]
lower.stt<-st.mod.t$conf.int[1]
stt.ci<-upper.stt-stt.est

# Direct estimation
exp.direct<-subset(exp,exp$direct==1|exp$direct==0)
direct.df<-summarySE(exp.direct, measurevar="direct", 
                     groupvars=c("ch1"),conf.interval = 0.90)

st.df<-data.frame(ch1=c("a","b"),N=c(NA,NA),direct=c(stc.est,stt.est),
                  sd=c(NA,NA),se=c(NA,NA),ci=c(stc.ci,stt.ci))
ict2.df<-rbind(direct.df,st.df)
ict2.df$mode<-c("Direct","Direct","List","List")

# plot2.2
f2.2.plot <- ggplot(ict2.df, 
                    aes(y=direct, x=ch1, group=mode, shape=mode))+
  scale_shape_manual(values =c(2,17))+
  scale_x_discrete(labels=c("COVID-19","Chinese Virus"))+
  theme_article()+geom_hline(yintercept = 0)
f2.2.plot <- f2.2.plot + geom_point(size=4, position=pos2) + 
  geom_errorbar(aes(ymin=direct-ci, ymax=direct+ci), width=0, colour="black", 
                position = pos2) + 
  labs(x=element_blank(),y="Estimated Proportion") + 
  theme(legend.title=element_blank(),
        axis.title.y=element_text(size=10),axis.text=element_text(size=10),
        legend.justification=c(1,1), legend.position=c(1,1))+
  guides(colour = FALSE, shape = guide_legend(reverse=T))
f2.2.plot

# combine plot 2.1 & plot 2.2
f2.plot<-ggarrange(f2.1.plot+theme(axis.text.x = element_blank(),
                                   axis.ticks.x = element_blank()),
                   f2.2.plot,
                   ncol = 1,nrow = 2)

f2.plot<-plot_grid(f2.plot)
f2.plot #Figure 2
# export Figure 2
ggsave(file = "f2.eps", width = 4.5, height = 3, dpi = 666)

####Figure 3####
##The “Chinese Virus” Label Treatment Effect

# DIM (bootstrap CI)
exp.covid<-subset(exp,ch==0)
exp.virus<-subset(exp,ch==1)

set.seed(123)
bootstrapped.c <- two.boot(exp.covid$ict.c[exp.covid$covid==1], 
                           exp.covid$ict.c[exp.covid$covid==0],
                           mean,
                           10000)
boot.mean_diff.covid <- data.frame(bootstrapped.c$t)
colnames(boot.mean_diff.covid) <- 'mean_diffs_c'

bootstrapped.t <- two.boot(exp.virus$ict.t[exp.virus$virus==1], 
                           exp.virus$ict.t[exp.virus$virus==0],
                           mean,
                           10000)

boot.mean_diff.virus <- data.frame(bootstrapped.t$t)
colnames(boot.mean_diff.virus) <- 'mean_diffs_v'
diffindiff <- boot.mean_diff.virus$mean_diffs_v-boot.mean_diff.covid$mean_diffs_c
mean.diffindiff <- mean(diffindiff)
lower.diffindiff <- quantile(diffindiff, probs = 0.05)
upper.diffindiff <- quantile(diffindiff, probs = 0.95)
cbind(lower.diffindiff, mean.diffindiff, upper.diffindiff)

# MLE
ml.results <- ictreg(ict~ch,
                     treat = "list",
                     data=exp,
                     J=3, method = "ml")

summary(ml.results)

exp.virus <- exp.covid <- exp
exp.virus[, "ch"] <- 1
exp.covid[, "ch"] <- 0

avg.pred.diff.mle <- predict(ml.results,
                             newdata = exp.virus, 
                             newdata.diff = exp.covid,
                             se.fit = TRUE, avg = TRUE, 
                             interval="confidence",
                             level=0.90)
avg.pred.diff.mle$fit
covid.mle.est<-avg.pred.diff.mle$fit[2,1]
virus.mle.est<-avg.pred.diff.mle$fit[1,1]
diff.mle.est<-avg.pred.diff.mle$fit[3,1]
covid.mle.upr<-avg.pred.diff.mle$fit[2,3]
covid.mle.lwr<-avg.pred.diff.mle$fit[2,2]
virus.mle.upr<-avg.pred.diff.mle$fit[1,3]
virus.mle.lwr<-avg.pred.diff.mle$fit[1,2]
diff.mle.upr<-avg.pred.diff.mle$fit[3,3]
diff.mle.lwr<-avg.pred.diff.mle$fit[3,2]

# NLS
nls.results <- ictreg(ict~ch,
                      treat = "list",
                      data=exp,
                      J=3, method = "nls")

summary(nls.results)

exp.virus <- exp.covid <- exp
exp.virus[, "ch"] <- 1
exp.covid[, "ch"] <- 0

avg.pred.diff.nls <- predict(nls.results,
                             newdata = exp.virus, 
                             newdata.diff = exp.covid,
                             se.fit = TRUE, avg = TRUE, 
                             interval="confidence",
                             level=0.90)

avg.pred.diff.nls$fit
covid.nls.est<-avg.pred.diff.nls$fit[2,1]
virus.nls.est<-avg.pred.diff.nls$fit[1,1]
diff.nls.est<-avg.pred.diff.nls$fit[3,1]
covid.nls.upr<-avg.pred.diff.nls$fit[2,3]
covid.nls.lwr<-avg.pred.diff.nls$fit[2,2]
virus.nls.upr<-avg.pred.diff.nls$fit[1,3]
virus.nls.lwr<-avg.pred.diff.nls$fit[1,2]
diff.nls.upr<-avg.pred.diff.nls$fit[3,3]
diff.nls.lwr<-avg.pred.diff.nls$fit[3,2]

## Combine bootmean, MLE, NLS plots [MAIN]
boot.mle.nls<-data.frame(est=c(stc.est,covid.mle.est,covid.nls.est,
                               stt.est,virus.mle.est,virus.nls.est,
                               mean.diffindiff,diff.mle.est,diff.nls.est),
                         upr=c(upper.stc,covid.mle.upr,covid.nls.upr,
                               upper.stt,virus.mle.upr,virus.nls.upr,
                               upper.diffindiff,diff.mle.upr,diff.nls.upr),
                         lwr=c(lower.stc,covid.mle.lwr,covid.nls.lwr,
                               lower.stt,virus.mle.lwr,virus.nls.lwr,
                               lower.diffindiff,diff.mle.lwr,diff.nls.lwr),
                         mtd=c("DIM","MLE","NLS","DIM","MLE","NLS","DIM","MLE","NLS"),
                         trt=c("COVID-19","COVID-19","COVID-19",
                               "Chinese Virus","Chinese Virus","Chinese Virus",
                               "Difference","Difference","Difference"))
boot.mle.nls$trt <- factor(boot.mle.nls$trt, levels=c("COVID-19","Chinese Virus","Difference"))

f3.plot <- ggplot(boot.mle.nls, 
                  aes(y=est, x=trt, group=mtd, shape=mtd))+
  scale_shape_manual(labels=c("DIM", "MLE", "NLS"),
                     values =myShape)+
  theme_article()+geom_hline(yintercept = 0)
f3.plot <- f3.plot + geom_point(size=4, position=pos3) + 
  geom_errorbar(limits, width=0, colour="black", position = pos3) + 
  labs(x=element_blank(),y="Estimated Proportion/Difference") + 
  theme(legend.title=element_blank(),
        axis.title.y=element_text(size=10),axis.text=element_text(size=10),
        legend.justification=c(0,0), legend.position=c(0,0))

f3.plot # Figure 3
# Export Figure 3
ggsave(file = "f3.eps", width = 6, height = 3, dpi = 666)

####Figure 4####
##Distribution Density of Blameworthiness
mean(exp$blame,na.rm=T)
sd(exp$blame,na.rm=T)
mean(exp$blame[exp$con==1],na.rm=T)
sd(exp$blame[exp$con==1],na.rm=T)
mean(exp$blame[exp$lib==1],na.rm=T)
sd(exp$blame[exp$lib==1],na.rm=T)

# Overall blame density
blamep<-ggplot(exp, aes(x=blame)) +
  geom_density(fill="gray32",color=NA)+
  labs(title ="Overall")+
  theme_article()
blamep<-blamep+theme(plot.title = element_text(size=10,hjust = 0.5))
blamep<-ggpar(blamep,font.tickslab=c(8,"plain","black"), 
              ylim = c(0.00, 0.02))
blamep

# by ideology
exp$ido1<-factor(exp$ido1)
is.factor(exp$ido1)
exp$ido2[exp$ido1=="1"]<-"Conservative"
exp$ido2[exp$ido1=="2"]<-"Moderate"
exp$ido2[exp$ido1=="3"]<-"Liberal"
exp$ido2[exp$ido1=="4"]<-"Moderate"

blamep.i<-ggplot(exp, aes(x=blame, fill=ido2)) +
  geom_density(alpha=0.7,color=NA)+
  scale_fill_manual(values=c("black","grey","white"))+
  guides(fill=guide_legend(title="Ideology"),color=FALSE)+
  labs(title ="Ideology")+
  theme_article()
blamep.i<-blamep.i+theme(plot.title = element_text(hjust = 0.5,size=10),
                         legend.title=element_blank(),
                         legend.text = element_text(size=8),
                         legend.justification=c(0.5,1),
                         legend.position=c(0.5,1),
                         legend.key = element_rect(colour = "black"))
blamep.i<-ggpar(blamep.i,font.tickslab=c(8,"plain","black"), 
                ylim = c(0.00, 0.02))
blamep.i

# combine overall and ideology subgroups
f4.plot<-ggarrange(blamep+theme(axis.title.y = element_blank(),
                                axis.title.x = element_blank()),
                   blamep.i+theme(axis.title.y = element_blank(),
                                  axis.title.x = element_blank(),
                                  axis.ticks.y = element_blank(),
                                  axis.text.y = element_blank()),
                   
                   ncol = 2,nrow = 1)
f4.plot<-plot_grid(f4.plot)
f4.plot #Figure 4

# export Figure 4
#options(bitmapType="cairo")
#pdf("f4.pdf", width = 6, height = 4)
#f4.plot
#dev.off()

####Figure 5####
##Heterogeneity in the “Chinese Virus” Label Treatment Effect

# Liberal
nls.lib <- ictreg(ict~ch,
                  treat = "list",
                  data=exp.lib,
                  J=3, method = "nls")

summary(nls.lib)

exp.virus.l <- exp.covid.l <- exp.lib
exp.virus.l[, "ch"] <- 1
exp.covid.l[, "ch"] <- 0

avg.pred.diff.lib <- predict(nls.lib,
                             newdata = exp.virus.l, 
                             newdata.diff = exp.covid.l,
                             se.fit = TRUE, avg = TRUE, 
                             interval="confidence",
                             level=0.90)

z=0.1955/0.1147
2*pnorm(-abs(z))

avg.pred.diff.lib$fit
avg.pred.diff.lib$se.fit
covid.lib.est<-avg.pred.diff.lib$fit[2,1]
virus.lib.est<-avg.pred.diff.lib$fit[1,1]
diff.lib.est<-avg.pred.diff.lib$fit[3,1]
covid.lib.upr<-avg.pred.diff.lib$fit[2,3]
covid.lib.lwr<-avg.pred.diff.lib$fit[2,2]
virus.lib.upr<-avg.pred.diff.lib$fit[1,3]
virus.lib.lwr<-avg.pred.diff.lib$fit[1,2]
diff.lib.upr<-avg.pred.diff.lib$fit[3,3]
diff.lib.lwr<-avg.pred.diff.lib$fit[3,2]

# Moderate
nls.mod <- ictreg(ict~ch,
                  treat = "list",
                  data=exp.mod,
                  J=3, method = "nls")

summary(nls.mod)

exp.virus.m <- exp.covid.m <- exp.mod
exp.virus.m[, "ch"] <- 1
exp.covid.m[, "ch"] <- 0

avg.pred.diff.mod <- predict(nls.mod,
                             newdata = exp.virus.m, 
                             newdata.diff = exp.covid.m,
                             se.fit = TRUE, avg = TRUE, 
                             interval="confidence",
                             level=0.90)

z=-0.0553/0.1752
2*pnorm(-abs(z))

avg.pred.diff.mod$fit
covid.mod.est<-avg.pred.diff.mod$fit[2,1]
virus.mod.est<-avg.pred.diff.mod$fit[1,1]
diff.mod.est<-avg.pred.diff.mod$fit[3,1]
covid.mod.upr<-avg.pred.diff.mod$fit[2,3]
covid.mod.lwr<-avg.pred.diff.mod$fit[2,2]
virus.mod.upr<-avg.pred.diff.mod$fit[1,3]
virus.mod.lwr<-avg.pred.diff.mod$fit[1,2]
diff.mod.upr<-avg.pred.diff.mod$fit[3,3]
diff.mod.lwr<-avg.pred.diff.mod$fit[3,2]

# Conservative
nls.con <- ictreg(ict~ch,
                  treat = "list",
                  data=exp.con,
                  J=3, method = "nls")

summary(nls.con)

exp.virus.c <- exp.covid.c <- exp.con
exp.virus.c[, "ch"] <- 1
exp.covid.c[, "ch"] <- 0

avg.pred.diff.con <- predict(nls.con,
                             newdata = exp.virus.c, 
                             newdata.diff = exp.covid.c,
                             se.fit = TRUE, avg = TRUE, 
                             interval="confidence",
                             level=0.90)

z=-0.0446/0.2103
2*pnorm(-abs(z))

avg.pred.diff.con$fit
covid.con.est<-avg.pred.diff.con$fit[2,1]
virus.con.est<-avg.pred.diff.con$fit[1,1]
diff.con.est<-avg.pred.diff.con$fit[3,1]
covid.con.upr<-avg.pred.diff.con$fit[2,3]
covid.con.lwr<-avg.pred.diff.con$fit[2,2]
virus.con.upr<-avg.pred.diff.con$fit[1,3]
virus.con.lwr<-avg.pred.diff.con$fit[1,2]
diff.con.upr<-avg.pred.diff.con$fit[3,3]
diff.con.lwr<-avg.pred.diff.con$fit[3,2]

## Combine ideology plots
ideo.df<-data.frame(est=c(covid.lib.est,virus.lib.est,diff.lib.est,
                          covid.mod.est,virus.mod.est,diff.mod.est,
                          covid.con.est,virus.con.est,diff.con.est),
                    upr=c(covid.lib.upr,virus.lib.upr,diff.lib.upr,
                          covid.mod.upr,virus.mod.upr,diff.mod.upr,
                          covid.con.upr,virus.con.upr,diff.con.upr),
                    lwr=c(covid.lib.lwr,virus.lib.lwr,diff.lib.lwr,
                          covid.mod.lwr,virus.mod.lwr,diff.mod.lwr,
                          covid.con.lwr,virus.con.lwr,diff.con.lwr),
                    mtd=c("a","b","c",
                          "a","b","c",
                          "a","b","c"),
                    trt=c("Liberal","Liberal","Liberal",
                          "Moderate","Moderate","Moderate",
                          "Conservative","Conservative","Conservative"))
ideo.df$trt <- factor(ideo.df$trt, 
                      levels=c("Liberal","Moderate","Conservative"))

f5.plot <- ggplot(ideo.df, 
                    aes(y=est, x=trt, group=mtd, shape=mtd))+
  scale_shape_manual(labels=c("COVID-19", "Chinese Virus", "Difference"),
                     values =myShape)+
  theme_article()+geom_hline(yintercept = 0)
f5.plot <- f5.plot + geom_point(size=4, position=pos3) + 
  geom_errorbar(limits, width=0, colour="black", position = pos3) + 
  labs(x=element_blank(),y="Proportion") + 
  theme(legend.title=element_blank(),axis.text=element_text(size=10),
        legend.position = "bottom")
f5.plot #Figure 5
# Export Figure 5
ggsave(file = "f5.eps", width = 6, height = 3, dpi = 666)

####Table 1####
##Item Count in the List Experiment

covid3<-table(exp$lc_c)
covid4<-table(exp$lt_c)
virus3<-table(exp$lc_t)
virus4<-table(exp$lt_t)

ict_stat<-matrix(NA,6,8)
colnames(ict_stat)<-c("Frequency","Proportion",
                      "Frequency","Proportion",
                      "Frequency","Proportion",
                      "Frequency","Proportion")
rownames(ict_stat)<-c("0","1","2","3","4","Sample size")
ict_stat[,1]<-c(covid3[1],covid3[2],covid3[3],covid3[4],NA,sum(covid3))
ict_stat[,2]<-c(covid3[1]/sum(covid3),covid3[2]/sum(covid3),
                covid3[3]/sum(covid3),covid3[4]/sum(covid3),NA,NA)
ict_stat[,3]<-c(covid4[1],covid4[2],covid4[3],covid4[4],covid4[5],sum(covid4))
ict_stat[,4]<-c(covid4[1]/sum(covid4),covid4[2]/sum(covid4),
                covid4[3]/sum(covid4),covid4[4]/sum(covid4),
                covid4[5]/sum(covid4),NA)
ict_stat[,5]<-c(virus3[1],virus3[2],virus3[3],virus3[4],NA,sum(virus3))
ict_stat[,6]<-c(virus3[1]/sum(virus3),virus3[2]/sum(virus3),
                virus3[3]/sum(virus3),virus3[4]/sum(virus3),NA,NA)
ict_stat[,7]<-c(virus4[1],virus4[2],virus4[3],virus4[4],virus4[5],sum(virus4))
ict_stat[,8]<-c(virus4[1]/sum(virus4),virus4[2]/sum(virus4),
                virus4[3]/sum(virus4),virus4[4]/sum(virus4),
                virus4[5]/sum(virus4),NA)

table1<-round(ict_stat,2)
print(xtable(x=table1,caption ="Item Count"))
write.table(table1, file = "table1.txt", sep = ",", quote = FALSE, row.names = T)

####Table 2####
##Estimation of Social Desirability Bias

Yc<-lm(Y.c~covid, data = exp)
summary(Yc)
Yt<-lm(Y.t~virus, data = exp)
summary(Yt)
Yc.l<-lm(Y.c~covid, data = exp.lib)
summary(Yc.l)
Yt.l<-lm(Y.t~virus, data = exp.lib)
summary(Yt.l)
Yc.c<-lm(Y.c~covid, data = exp.con)
summary(Yc.c)
Yt.c<-lm(Y.t~virus, data = exp.con)
summary(Yt.c)
Yc.m<-lm(Y.c~covid, data = exp.mod)
summary(Yc.m)
Yt.m<-lm(Y.t~virus, data = exp.mod)
summary(Yt.m)

# Overall
covid.coef<-summary(Yc)$coefficients[2,1]
covid.sd<-summary(Yc)$coefficients[2,2]
covid.p<-summary(Yc)$coefficients[2,4]
virus.coef<-summary(Yt)$coefficients[2,1]
virus.sd<-summary(Yt)$coefficients[2,2]
virus.p<-summary(Yt)$coefficients[2,4]
# Lib
covid.coef.l<-summary(Yc.l)$coefficients[2,1]
covid.sd.l<-summary(Yc.l)$coefficients[2,2]
covid.p.l<-summary(Yc.l)$coefficients[2,4]
virus.coef.l<-summary(Yt.l)$coefficients[2,1]
virus.sd.l<-summary(Yt.l)$coefficients[2,2]
virus.p.l<-summary(Yt.l)$coefficients[2,4]
# Con
covid.coef.c<-summary(Yc.c)$coefficients[2,1]
covid.sd.c<-summary(Yc.c)$coefficients[2,2]
covid.p.c<-summary(Yc.c)$coefficients[2,4]
virus.coef.c<-summary(Yt.c)$coefficients[2,1]
virus.sd.c<-summary(Yt.c)$coefficients[2,2]
virus.p.c<-summary(Yt.c)$coefficients[2,4]
# Mod
covid.coef.m<-summary(Yc.m)$coefficients[2,1]
covid.sd.m<-summary(Yc.m)$coefficients[2,2]
covid.p.m<-summary(Yc.m)$coefficients[2,4]
virus.coef.m<-summary(Yt.m)$coefficients[2,1]
virus.sd.m<-summary(Yt.m)$coefficients[2,2]
virus.p.m<-summary(Yt.m)$coefficients[2,4]

table2<-matrix(NA,8,6)
colnames(table2)<-c("Overt Perceived Threat",
                    "Social Desirability Effect",
                    "p-value","Overt Perceived Threat",
                    "Social Desirability Effect",
                    "p-value")
rownames(table2)<-c("Overall","S.E.",
                    "Liberal","S.E.","Conservative","S.E.",
                    "Moderate","S.E.")
table2[1,]<-c(mean(exp$direct_c,na.rm=T),covid.coef,covid.p,
              mean(exp$direct_t,na.rm=T),virus.coef,virus.p)
table2[2,]<-c(sd(exp$direct_c,na.rm=T),covid.sd,NA,
              sd(exp$direct_t,na.rm=T),virus.sd,NA)
table2[3,]<-c(mean(exp.lib$direct_c,na.rm=T),covid.coef.l,covid.p.l,
              mean(exp.lib$direct_t,na.rm=T),virus.coef.l,virus.p.l)
table2[4,]<-c(sd(exp.lib$direct_c,na.rm=T),covid.sd.l,NA,
              sd(exp.lib$direct_t,na.rm=T),virus.sd.l,NA)
table2[5,]<-c(mean(exp.con$direct_c,na.rm=T),covid.coef.c,covid.p.c,
              mean(exp.con$direct_t,na.rm=T),virus.coef.c,virus.p.c)
table2[6,]<-c(sd(exp.con$direct_c,na.rm=T),covid.sd.c,NA,
              sd(exp.con$direct_t,na.rm=T),virus.sd.c,NA)
table2[7,]<-c(mean(exp.mod$direct_c,na.rm=T),covid.coef.m,covid.p.m,
              mean(exp.mod$direct_t,na.rm=T),virus.coef.m,virus.p.m)
table2[8,]<-c(sd(exp.mod$direct_c,na.rm=T),covid.sd.m,NA,
              sd(exp.mod$direct_t,na.rm=T),virus.sd.m,NA)

table2<-round(table2,2)
print(xtable(x=table2))
write.table(table2, file = "table2.txt", sep = ",", quote = FALSE, row.names = T)

####Table 3####
##The Federal Government’s Blameworthiness

blame<-lm(blame~ch,data=exp)
summary(blame)
blame.l<-lm(blame~ch,data=exp.lib)
summary(blame.l)
blame.c<-lm(blame~ch,data=exp.con)
summary(blame.c)
blame.m<-lm(blame~ch,data=exp.mod)
summary(blame.m)

# Export Table 3
stargazer(blame,blame.l,blame.c,blame.m,
          style = "apsr",type="html",out="table3.htm",
          covariate.labels=c("Chinese Virus"),
          column.labels=c("Overall",
                          "Liberal","Conservative","Moderate"),
          dep.var.labels.include = FALSE,
          intercept.bottom=TRUE,
          keep.stat = c("n","rsq"),
          align=TRUE,no.space=TRUE)


